Xe & XeF2 notes and plotting

20/04/21

For manuscript notes & figures.

Based on http://127.0.0.1:8888/lab/tree/Scratch/XeF2-notebooks/xe-xef2_comp_120221-dist.ipynb and SO stuff from http://127.0.0.1:8888/lab/tree/Scratch/XeF2-notebooks/XeF2_ePS-expt_comp_121120_4d-dev.ipynb

  • Currently running on Stimpy
  • Use iPyPublish for final output (bemo only at the momemnt). Outputting selected cells & figures.

CURRENT VERSION (21/04/21): copied in some Lyx LaTex, set cell metadata and tested. Working... not very slick. Need to rerun for PDF figures.


12/02/21

Comparison of Xe & XeF2 data + experiments. No SO treatment herein.

Basic conclusions: atomic Xe parameters (cross-sections and $\beta$) look broadly similar to XeF2 results, except:

  • Much broader feature in the XS.
  • No additional structures in the $\beta$s, i.e. it's much smoother, for the atomic case. This confirms that the structures in the molecular case are 'real' molecular features (at least as far as the fixed-nuclei, single-electron calculations go), but whether they persist in realtiy and are experimentally observable is a good question.

For XeF2 only plus SO treatment, see previous notebook XeF2_ePS-expt_comp_271020_4d_v111120-dist.html.

Intro notes

Setup

In [1]:
!hostname
Stimpy
In [2]:
!conda env list
# conda environments:
#
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\dataVis3D
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\dataVis3D-yt
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSpkgTest2
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\fibre-sim
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ipykernel_py2
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\pkgTest
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\pypi-test
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\seabreeze-bk
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\test
base                     C:\ProgramData\Anaconda3
ePSdev                   C:\Users\femtolab\.conda\envs\ePSdev
epsInstallTest           C:\Users\femtolab\.conda\envs\epsInstallTest
epsdev                *  C:\Users\femtolab\.conda\envs\epsdev
epsman                   C:\Users\femtolab\.conda\envs\epsman
seabreeze                C:\Users\femtolab\AppData\Local\conda\conda\envs\seabreeze
tensorflow               C:\Users\femtolab\AppData\Local\conda\conda\envs\tensorflow

In [3]:
import sys
import os
from pathlib import Path
import numpy as np
# import epsproc as ep
import xarray as xr

import matplotlib.pyplot as plt

from datetime import datetime as dt
timeString = dt.now()
In [4]:
# For module testing, include path to module here, otherwise use global installation
if sys.platform == "win32":
    modPath = r'D:\code\github\ePSproc'  # Win test machine
    winFlag = True
else:
    modPath = r'/home/femtolab/github/ePSproc/'  # Linux test machine
    winFlag = False
    
sys.path.append(modPath)
import epsproc as ep

# Plotters
from epsproc.plot import hvPlotters

# Multijob class dev code
from epsproc.classes.multiJob import ePSmultiJob
* pyevtk not found, VTK export not available. 
In [5]:
hvPlotters.setPlotters(width = 1000, snsStyle="whitegrid")
# import bokeh
# import holoviews as hv
# hv.extension('bokeh')


# For class, above settings don't take, not sure why, something to do with namespaces/calling sequence?
# Overriding snsStyle does work however... although NOT CONSISTENTLY????
# AH, looks like ordering matters - set_style LAST (.set seems to override)
import seaborn as sns

sns.set(rc={'figure.figsize':(10,6)})  # Set figure size in inches
sns.set_context("paper")
sns.set_style("whitegrid")  # Set plot style
sns.set_palette("Paired")   # Set colour mapping

# sns.set_theme(context='paper', style='whitegrid')  # Not valid for v0.9.0
In [6]:
# Try direct fig type setting for PDF output figs
from IPython.display import set_matplotlib_formats
# set_matplotlib_formats('png', 'pdf')
set_matplotlib_formats('svg', 'pdf')

Load data

In [7]:
# # Scan for subdirs, based on existing routine in getFiles()

# fileBase = r'D:\projects\ePolyScat\xef2\XeF2_highRes_wf'  # Data dir on Stimpy
# fileBase = r'D:\VMs\Share\ePSshare\xe\Xe_wf'  # Data dir on Stimpy
# jobDirs = [r'D:\VMs\Share\ePSshare\xe\Xe_wf\orb19_4d']  # Check OK for single dir specified method

# Set explicit list of dirs
jobDirs = [r'D:\VMs\Share\ePSshare\xe\Xe_wf\orb19_4d', r'D:\projects\ePolyScat\xef2\XeF2_highRes_wf\orb21_A1G',
          r'D:\projects\ePolyScat\xef2\XeF2_highRes_wf\orb22_E1G', r'D:\projects\ePolyScat\xef2\XeF2_highRes_wf\orb24_E2G']
In [8]:
data = ePSmultiJob(jobDirs = jobDirs, verbose = 0)
In [9]:
# keys = [4,5,6]  # Set for 4d data only
# data.scanFiles(keys = keys)

data.scanFiles()  # All files
data.jobsSummary()
*** Warning: Missing records, expected 80, found 60.
*** Warning: Found 20 blank sets of matrix elements, symmetries ['HU']
*** Warning: Missing records, expected 80, found 60.
*** Warning: Found 20 blank sets of matrix elements, symmetries ['HU']
*** Warning: Missing records, expected 80, found 60.
*** Warning: Found 20 blank sets of matrix elements, symmetries ['HU']
*** Warning: Missing records, expected 80, found 60.
*** Warning: Found 20 blank sets of matrix elements, symmetries ['HU']
*** Warning: Missing records, expected 80, found 60.
*** Warning: Found 20 blank sets of matrix elements, symmetries ['HU']
*** Warning: Missing records, expected 80, found 60.
*** Warning: Found 20 blank sets of matrix elements, symmetries ['HU']
*** Warning: Missing records, expected 80, found 60.
*** Warning: Found 20 blank sets of matrix elements, symmetries ['HU']
*** Warning: Missing records, expected 80, found 60.
*** Warning: Found 20 blank sets of matrix elements, symmetries ['HU']
*** Warning: Missing records, expected 80, found 60.
*** Warning: Found 20 blank sets of matrix elements, symmetries ['HU']
*** Warning: Missing records, expected 80, found 60.
*** Warning: Found 20 blank sets of matrix elements, symmetries ['HU']
*** Warning: Missing records, expected 80, found 60.
*** Warning: Found 20 blank sets of matrix elements, symmetries ['HU']
*** Warning: Missing records, expected 80, found 60.
*** Warning: Found 20 blank sets of matrix elements, symmetries ['HU']
*** Warning: Missing records, expected 80, found 60.
*** Warning: Found 20 blank sets of matrix elements, symmetries ['HU']
*** Warning: Missing records, expected 80, found 60.
*** Warning: Found 20 blank sets of matrix elements, symmetries ['HU']
*** Warning: Missing records, expected 80, found 60.
*** Warning: Found 20 blank sets of matrix elements, symmetries ['HU']
Found 4 directories, with 95 files.

*** Job orb19_4d details
Key: orb19_4d
Dir D:\VMs\Share\ePSshare\xe\Xe_wf\orb19_4d, 15 file(s).
{   'batch': 'ePS xe, batch Xe_wf, orbital orb19_4d',
    'event': 'orb 19 ionization (4d, S/Ag), wavefunction run.',
    'orbE': -72.52106488190199,
    'orbLabel': '4d, S/Ag'}

*** Job orb21_A1G details
Key: orb21_A1G
Dir D:\projects\ePolyScat\xef2\XeF2_highRes_wf\orb21_A1G, 30 file(s).
{   'batch': 'ePS XeF2_2020, batch XeF2_highRes_wf, orbital orb21_A1G',
    'event': 'orb 21 ionization (Xe 4d, A1G/SG), sph grid. Inputs based on '
             'original 2019 calcs, now with chunking for higher E resolution.',
    'orbE': -76.581003676086,
    'orbLabel': 'Xe 4d, A1G/SG'}

*** Job orb22_E1G details
Key: orb22_E1G
Dir D:\projects\ePolyScat\xef2\XeF2_highRes_wf\orb22_E1G, 30 file(s).
{   'batch': 'ePS XeF2_2020, batch XeF2_highRes_wf, orbital orb22_E1G',
    'event': 'orb 22/23 ionization (Xe 4d, E1G/PG), sph grid. Inputs based on '
             'original 2019 calcs, now with chunking for higher E resolution.',
    'orbE': -76.27895729126399,
    'orbLabel': 'Xe 4d, E1G/PG'}

*** Job orb24_E2G details
Key: orb24_E2G
Dir D:\projects\ePolyScat\xef2\XeF2_highRes_wf\orb24_E2G, 20 file(s).
{   'batch': 'ePS XeF2_2020, batch XeF2_highRes_wf, orbital orb24_E2G',
    'event': 'orb 24/25 ionization (Xe 4d, E2G/DG), sph grid. Inputs based on '
             'original 2019 calcs, now with chunking for higher E resolution.',
    'orbE': -75.67486452162,
    'orbLabel': 'Xe 4d, E2G/DG'}

(Note: for the atomic case, there is an HU continuum symmetry allowed, but the matrix elements are, effectively, zero - this triggers the warning above.)

In [10]:
# Change labels for plotting
# data.data['orb19_4d']['jobNotes']['orbLabel'] = 'Xe atomic 4d, S/Ag'

for key in data.data.keys():
    data.data[key]['jobNotes']['orbLabel'] = data.data[key]['jobNotes']['orbLabel'].split(',')[-1]
#     print(data.data[key]['jobNotes']['orbLabel'].split(',')[-1])
In [11]:
# Set branching ratios

# Stack XS data to new data structure
# NOTE: this is not quite correct, since it forces all data to same Ehv axis, but OK for a quick hack, and easy to pull out branching ratios
dsXS = xr.Dataset()  # Set blank dataset, this is easier for stacking, probably
lText = []

for key in data.data.keys():
    if key is not 'orb19_4d':  # Skip Xe atomic case
        subset = data.data[key]['XS'].sel({'Sym':'All', 'XC':'SIGMA'})  # Set XS data, all syms only

        # NOTE currently missing full dataset resolution for orb24, so try interp. (2eV not 1eV step size)
        # Note dropna to ensure no NaNs, see http://xarray.pydata.org/en/latest/interpolation.html#interpolating-arrays-with-nan
        if key == 'orb24_E2G':
            subset = subset.dropna('Eke').interp(Eke = dsXS.Eke, method = 'cubic')

    #     dsXS[key] = subset.copy()  # Set .copy() for safety here!
        dsXS[data.data[key]['jobNotes']['orbLabel']]  = subset.copy()  # Set .copy() for safety here!


# Convert to Xarray & normalise
dsXS = dsXS.to_array().rename({'variable':'channel'}).squeeze()
dsXS['total'] = dsXS.sum('channel')  # Sum over channels

# Normalise...
dsXS = dsXS/dsXS['total']

Experimental data

  • Set 1: RF processed results, 02/10/20
  • Set 2: DH ion yields, 27/10/20
In [12]:
# Load experimental data
dataPathExpt = Path(r'D:\projects\XeF2_Soleil_2019\RF_data_analysis_021020')

# Set data attribs in dict similar to ePS results structure
dataFiles = {}
dataFiles['SIGMA'] = {'Data':'XS', 'File':r'XeF2_Xe4d_4d32_4d52_cross_sec_all_photon_energies_02102020.dat', 
                      'States':['$4d_{3/2}$', '$4d_{5/2}$']}
dataFiles['BETA'] = {'Data':'Beta', 'File':r'XeF2_Xe4d_beta_all_photon_energies_02102020.dat', 
                     'States':['$\Pi_{1/2} (4d_{3/2})$', '$\Delta_{3/2} (4d_{3/2})$', '$\Sigma_{1/2} (4d_{5/2})$', 
                               '$\Pi_{3/2} (4d_{5/2})$', '$\Delta_{5/2} (4d_{5/2})$']}

dataFiles['BR-All'] = {'Data':'Branching ratios', 'File':r'XeF2_Xe4d_branching_all_photon_energies_all_states_01112020.dat', 
                      'States':['$\Pi_{1/2} (4d_{3/2})$', '$\Delta_{3/2} (4d_{3/2})$', '$\Sigma_{1/2} (4d_{5/2})$', 
                               '$\Pi_{3/2} (4d_{5/2})$', '$\Delta_{5/2} (4d_{5/2})$']}
dataFiles['BR-SO-summed'] = {'Data':'Branching ratios SO summed', 'File':r'XeF2_Xe4d_branching_all_photon_energies_SO_av_01112020.dat', 
                     'States':['$\Pi$', '$\Delta$', '$\Sigma$']}

# Update with ion data
dataFiles2 = {}
dataFiles2['ION-low'] = {'Data':'XS', 'File':r'XeF2_ion_yield_low_energy_cal.txt'}
dataFiles2['ION-high'] = {'Data':'XS', 'File':r'XeF2_ion_yield_high_energy_cal.txt'}
                         
# Update with branching ratios
# dataFilesBR = {}
# dataFilesBR['All'] = {'Data':'Branching ratios', 'File':r'XeF2_Xe4d_branching_all_photon_energies_all_states_01112020.dat', 
#                       'States':['\pi1/2 (4d3/2)', '\delta3/2 (4d3/2)', '\sigma1/2 (4d5/2)', '\pi3/2 (4d5/2)', '\delta5/2 (4d5/2)']}
# dataFiles['SO-summed'] = {'Data':'Branching ratios SO summed', 'File':r'XeF2_Xe4d_branching_all_photon_energies_SO_av_01112020.dat', 
#                      'States':['\pi', '\delta', '\sigma']}
In [13]:
# Read data files and convert to Xarray
# 27/10/20 added quick hack to set 2nd array for ion yield data
dataList = []
dataList2 = []
for key in dataFiles:
    dataIn = np.loadtxt(dataPathExpt/dataFiles[key]['File'])
    
    # Convert to Xarray
    dataXr = xr.DataArray(dataIn[:,1:], dims=('Ehv','State'), coords={'Ehv':dataIn[:,0], 'State':dataFiles[key]['States'][0:dataIn.shape[1]-1]})
    dataXr.attrs = dataFiles[key]
    dataXr.attrs['dataPath'] = dataPathExpt
    dataXr.name = key
    dataList.append(dataXr)

# Stack to Xarray
dataExpt = xr.concat(dataList, "XC")
dataExpt['XC'] = list(dataFiles.keys())

for key in dataFiles2:
    dataIn = np.loadtxt(dataPathExpt/dataFiles2[key]['File'])
    
    # Convert to Xarray
    dataXr = xr.DataArray(dataIn[:,1].squeeze(), dims=('Ehv'), coords={'Ehv':dataIn[:,0]})  # Only 1D datasets in this case
    dataXr.attrs = dataFiles2[key]
    dataXr.attrs['dataPath'] = dataPathExpt
    dataXr.name = key
    dataList2.append(dataXr)

dataExpt2 = xr.concat(dataList2, "XC")
dataExpt2['XC'] = list(dataFiles2.keys())

Comparison plots...

Here set Seaborn pallete for colour-mapping, and include labels etc. - aim to be publication ready.

In [15]:
# sns.color_palette("cubehelix")
# sns.set_palette("Reds")
# sns.set_palette("Paired")

# Compare with computational results
marker = None #'x'
Etype = 'Ehv'
Erange = [70, 255]

pType = 'SIGMA'
pltObj, lTextComp = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)
lText = lTextComp.copy()

# Add expt data - cross-secions
scale = dataExpt.sel({'XC':pType}).max() / data.dataSets['orb19_4d']['XS'].sel({'XC':pType, 'Type':'L'}).max()  # Set scaling to match computational data
(dataExpt.sel({'XC':pType})/scale).dropna('State').plot.line(x='Ehv', marker=marker, ls='--');

# Add expt data - ion yields
ionData = 'ION-high'
scale = dataExpt2.sel({'XC':ionData}).max() / data.dataSets['orb21_A1G']['XS'].sel({'XC':pType, 'Type':'L'}).max()  # Set scaling to match computational data
(dataExpt2.sel({'XC':ionData})/scale).plot.line(x='Ehv', marker=marker, ls='--');

# Manual legend fix
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
lText.extend([ionData])
plt.legend(lText)

# Fix axis labels
if pType == 'SIGMA':
    plt.ylabel('Cross-section/Mb')
#     plt.title('Cross-sections (normalised to computational results)')
    plt.title('')
else:
    plt.ylabel(r"$\beta_{LM}$")
#     plt.title(r"$\beta_{LM}$")
    
# Fix plot x-axis
plt.xlim(Erange);

# plt.yscale("log")  # Use log scale y-axis?
2021-04-22T09:40:56.785008 image/svg+xml Matplotlib v3.3.1, https://matplotlib.org/

Comparison of experimental (dashed lines) and computational (solid lines) cross-sections. Computational cross-sections are in Mb, and experimental results are scaled to match the peak value.

In [16]:
# Beta comparison plot over orbs
pType = 'BETA'
pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)

# Add expt data
dataExpt.sel({'XC':pType}).dropna('State').plot.line(x='Ehv', marker=marker, ls='--');

# Manual legend fix
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText, loc='upper right')

# Fix axis labels
if pType == 'SIGMA':
    plt.ylabel('XS/Mb')
    plt.title('XS (normalised to computational results)')
else:
    plt.ylabel(r"$\beta_{LM}$")
#     plt.title(r"$\beta_{LM}$")
    plt.title('')
    
# Fix plot x-axis
plt.xlim(Erange);
2021-04-22T09:41:22.903583 image/svg+xml Matplotlib v3.3.1, https://matplotlib.org/

Comparison of experimental (dashed lines) and computational (solid lines) $\beta$ parameters.

In [17]:
# Beta comparison plot over orbs
Erange = [120, 200]

pType = 'BETA'
pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)

# Add expt data
dataExpt.sel({'XC':pType}).dropna('State').plot.line(x='Ehv', marker=marker, ls='--');

# Manual legend fix
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText, loc='upper right')

# Fix axis labels
if pType == 'SIGMA':
    plt.ylabel('XS/Mb')
    plt.title('XS (normalised to computational results)')
else:
    plt.ylabel(r"$\beta_{LM}$")
    plt.title(r"$\beta_{LM}$")
    
# Fix plot x-axis
plt.xlim(Erange);
2021-04-22T09:41:23.461708 image/svg+xml Matplotlib v3.3.1, https://matplotlib.org/
In [18]:
# Compare with computational results
Erange = [75, 250]
pType = 'BR-SO-summed'
marker = 'x'

# pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)

dsXS.swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L').plot.line(x='Ehv');

# Add expt data - cross-secions
# scale = dataExpt.sel({'XC':pType}).max() / data.dataSets['orb21']['XS'].sel({'XC':pType, 'Type':'L'}).max()  # Set scaling to match computational data
scale = 1
(dataExpt.sel({'XC':pType})/scale).dropna('State').plot.line(x='Ehv', marker=marker, ls=':');

# Manual legend fix
lText = lTextComp.copy()[1:]
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data);
plt.legend(lText, loc='upper right');

plt.ylabel("Branching ratio")
plt.title("");
2021-04-22T09:41:23.958487 image/svg+xml Matplotlib v3.3.1, https://matplotlib.org/

Comparison of experimental (dashed lines) and computational (solid lines) branching ratios. Experimental results are summed over spin-orbit bands to allow direct comparison with the computational results.

Spin-orbit

In [19]:
from epsproc.geomFunc.geomUtils import genllLList
from epsproc.geomFunc.geomCalc import w3jTable

# Gen QNs for specific (L,S) case
L = 2
S = 0.5

Llist = np.array([[L,L+S,S], [L, L, S], [L,L-S,S]])   # Note this needs to be 2D array in current form of function
QNs = genllLList(Llist, uniqueFlag = False)

# Then calc 3js....
backend = 'sympy'
form = 'xdaLM'  # '2d'  # 'xdaLM'  # xds
nonzeroFlag = True

dlist = ['L', 'J', 'S', 'Lambda', 'Omega', 'Sigma']  # Set dims for reference

thrj = w3jTable(QNs = QNs, nonzeroFlag = nonzeroFlag, form = form, dlist = dlist, backend = backend)

# And primed terms (will be identical at this point, but set dims for multiplication later)
thrjP = w3jTable(QNs = QNs, nonzeroFlag = nonzeroFlag, form = form, dlist = [item + 'p' for item in dlist], backend = backend)
In [20]:
# Reformate by comsolidating +/- terms as modulation for XS and square
thrjUS = thrj.unstack('mSet').fillna(0)
# thrjXSmod = 2*(thrjUS.sum('Sigma').where((thrjUS['Lambda']>-0.5)&(thrjUS['Omega']<0.5), drop=True))**2  # Note Lambda & Omega anti-phase!
thrjXSmod = 2*(thrjUS.sum('Sigma').where((thrjUS['Lambda']>-0.5)&(thrjUS['Omega']<0.5), drop=True))**2  # Note Lambda & Omega anti-phase!
thrjXSmod['Omega'] = np.abs(thrjXSmod['Omega'])

pdTable, _ = ep.multiDimXrToPD(thrjXSmod, colDims = 'Lambda')
pdTable
Out[20]:
Lambda 0.0 1.0 2.0
J L Omega S
1.5 2.0 2.5 0.5 0.0 0.000000 0.000000
1.5 0.5 0.0 0.100000 0.400000
0.5 0.5 0.2 0.300000 0.000000
2.5 2.0 2.5 0.5 0.0 0.000000 0.333333
1.5 0.5 0.0 0.266667 0.066667
0.5 0.5 0.2 0.133333 0.000000
In [21]:
# SO branching ratios

# 11/11/20 - converted to function

def simulateSOBR(dsXS, thrjSOTerm, BRtype = 'All', coupling = 'Lambda', stateLabels = ['$\Sigma$', '$\Pi$', '$\Delta$']):
    # Set Lambda terms
    # dsXS['Lambda'] = dsXS['channel']
    # dsXS['Lambda'] = [0.0, 1.0, 2.0]
    
    dsXSO = dsXS.assign_coords({coupling:('channel',thrjSOTerm[coupling])}).swap_dims({'channel':coupling})  # Assign Lambda for multiplication
    
    dsXSO = dsXSO.assign_coords({f'{coupling}-Labels':(coupling, np.asarray(stateLabels)[dsXSO[coupling].astype(int).data])})  # Assign labels, based on int values as indexers (CRUDE!)


        
    dsXSO = dsXSO * thrjSOTerm 
    # dsXS

    # Convert to branching ratios (total)
    if BRtype == 'All':
#         dsXSO['total'] = dsXSO.sum(['Lambda', 'Omega', 'lSet'])  # Sum over channels
        dsXSO['total'] = dsXSO.sum(thrjSOTerm.dims)
    
    elif BRtype == 'J':
        dsXSO['total'] = dsXSO.sum(['Lambda', 'Omega'])  # Sum over channels, except J - GIVES same result? 
        
    elif BRtype == 'None':
        dsXSO['total'] = 1.0
        
    else:
        print(f'BRtype={BRtype} not recognised.')
        return

    # # Set branching ratios
    dsXSO = dsXSO/dsXSO['total']

    return dsXSO

State selected Assuming both $\Lambda$ and $\Omega$ are defined... as per experimental assignments (?)

Label levels by ($\Omega$, $\Lambda$), with no additional summation over terms. Q: is this justified/realistic?

In [22]:
pType = 'BR-All'

# Selected term(s) only - LOOKS GOOD FOR coupling='Omega' case, but opposite phases for coupling='Lambda'
dsXSO = simulateSOBR(dsXS, thrjXSmod, BRtype = 'None', coupling='Omega')  # This looks good
# dsXSO = simulateSOBR(dsXS, np.abs(thrjXSmodCoherent), BRtype = 'None', coupling='Lambda')  # Testing coherent case

# Unnorm
# dsXSO *= dsXSO['total']

# Selected term(s) only
J=1.5
OLset = [[2,1.5], [1,0.5]] # [Omega, Lambda] pairs to match expt.
# OLset = [[1,0.5],[2,1.5]] # [Omega, Lambda] pairs to match expt.
# OLset = [[1.5,2], [0.5,1]]  # Check ordering OK! 

# J=2.5
# OLset = [[2,2.5], [1,1.5], [0,0.5]] # [Omega, Lambda] pairs to match expt.

# subset = xr.DataArray()
subset = []
for OL in OLset:
    subset.append((dsXSO).swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J, Lambda=OL[0], Omega=OL[1]))

subsetXS = xr.concat(subset, dim='Omega').squeeze()

# Renorm
subsetXS = subsetXS/subsetXS.sum('Omega')

subsetXS.plot.line(x=Etype)

# Selected states only
if J == 1.5:
    statesPlot = ['$\\Delta_{3/2} (4d_{3/2})$', '$\\Pi_{1/2} (4d_{3/2})$']
else:
    statesPlot = ['$\\Delta_{5/2} (4d_{5/2})$', '$\\Pi_{3/2} (4d_{5/2})$', '$\\Sigma_{1/2} (4d_{5/2})$']
    
(dataExpt.sel({'XC':pType, 'State':statesPlot})/scale).plot.line(x='Ehv', marker=marker, ls=':');


# Manual legend fix
# lText = list(subsetXS['Omega'].data)
# lText = list(OLset)
lText = [f"{item.split(' ')[0]}$(sim)" for item in statesPlot]
# lText = list(subsetXS['Lambda'].data)
lText.extend(statesPlot)
# lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText, loc='lower right')

plt.ylabel("Branching ratio")
# plt.title(f'Simulated branching ratios, J={J}, Hund\'s case (c), with subselection (no summation)')
plt.title("");

plt.ylim([0.3, 0.68]);  # Set ylims to avoid large peaks
2021-04-22T09:41:24.502830 image/svg+xml Matplotlib v3.3.1, https://matplotlib.org/

Comparison of experimental (dashed lines) and computational (solid lines) branching ratios for $J=3/2$.

In [23]:
# Selected term(s) only
# Selected term(s) only - LOOKS GOOD FOR coupling='Omega' case, but opposite phases for coupling='Lambda'
# For J = 5/2 'Lambda' case is possibly better, at least for state ordering? DO NOT UNDERSTAND THIS YET!
# scale = 1
dsXSO = simulateSOBR(dsXS, thrjXSmod, BRtype = 'None', coupling='Omega')

# Unnorm
# dsXSO *= dsXSO['total']

# Selected term(s) only
# J=1.5
# OLset = [[2,1.5], [1,0.5]] # [Omega, Lambda] pairs to match expt.

J=2.5
# OLset = [[2,2.5], [1,1.5], [0,0.5]] # [Omega, Lambda] pairs to match expt.
OLset = [[0,0.5], [1,1.5], [2,2.5]] # [Omega, Lambda] pairs to match expt - REVERSED STATE ORDERING, justified by ion state energy ordering?

# subset = xr.DataArray()
subset = []
for OL in OLset:
    subset.append((dsXSO).swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J, Lambda=OL[0], Omega=OL[1]))

subsetXS = xr.concat(subset, dim='Omega').squeeze()

# Renorm
subsetXS = subsetXS/subsetXS.sum('Omega')

subsetXS.plot.line(x=Etype)

# Selected states only
if J == 1.5:
    statesPlot = ['$\\Delta_{3/2} (4d_{3/2})$', '$\\Pi_{1/2} (4d_{3/2})$']
else:
    statesPlot = ['$\\Delta_{5/2} (4d_{5/2})$', '$\\Pi_{3/2} (4d_{5/2})$', '$\\Sigma_{1/2} (4d_{5/2})$']
    
(dataExpt.sel({'XC':pType, 'State':statesPlot})/scale).plot.line(x='Ehv', marker=marker, ls=':');


# Manual legend fix
# lText = list(subsetXS['Omega'].data)
# lText = list(OLset)
lText = [f"{item.split(' ')[0]}$(sim)" for item in statesPlot]
# lText = list(subsetXS['Lambda'].data)
lText.extend(statesPlot)
# lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
# plt.legend(lText, loc='upper right')
plt.legend(lText, loc='lower right')

plt.ylabel("Branching ratio")
# plt.title(f'Simulated branching ratios, J={J}, Hund\'s case (c), with subselection (no summation)')
plt.title('');

plt.ylim([0.1, 0.5]);  # Set ylims to avoid large peaks
2021-04-22T09:41:24.996944 image/svg+xml Matplotlib v3.3.1, https://matplotlib.org/

Comparison of experimental (dashed lines) and computational (solid lines) branching ratios for $J=5/2$.

Full info dump

XS and $\beta$ per orbital & symmetry

Interactive^ plots per orbital & symmetry. In these plots the solid lines show the 'mixed' gauge results, and dashed lines show length and velocity gauges as bounds on the values.

^ Use controls in top right of plot. Also data sets can be turned on/off using the legend entries.

Note that Xe atomic case is the first plot below, ignore the out-of-bounds $\beta$s for (null) HU symmetry case, these are artefacts. Also of note in this case is that the individual symmetry channels show very little variation with energy, hence most of the overall (symmetry-summed) structure is form intra-channel effects (interferences).

In [24]:
# XC & Betas
data.plotGetCro(pType='SIGMA', Etype = Etype, Erange = Erange, backend = 'hv') 
*** 4d, S/Ag
*** Xe 4d, A1G/SG
*** Xe 4d, E1G/PG
*** Xe 4d, E2G/DG

Matrix elements

As previously for XeF2 case; for atomic Xe there are 5 degenerate channels (4d components) labelled by 'it'.

In [25]:
data.lmPlot(dataType = 'matE', Erange = Erange, Etype = Etype, thres = 0.1, logFlag = True, selDims = {'Type':'L', 'it':1}, fillna = True, cmap = "vlag", figsize = (15,8))   

# cmap = "cubehelix" )  # Quite good for showing resonance features.
Plotting data Xe_wf.orb19_4d_E1.0_15.0_286.0eV.inp.out, pType=a, thres=0.1, with Seaborn
Plotting data XeF2_highRes_wf.orb21_A1G_E1.0_10.0_141.0eV.inp.out, pType=a, thres=0.1, with Seaborn
Plotting data XeF2_highRes_wf.orb22_E1G_E1.0_10.0_141.0eV.inp.out, pType=a, thres=0.1, with Seaborn
Plotting data XeF2_highRes_wf.orb24_E2G_E1.0_10.0_141.0eV.inp.out, pType=a, thres=0.1, with Seaborn
2021-04-22T09:41:36.751070 image/svg+xml Matplotlib v3.3.1, https://matplotlib.org/
2021-04-22T09:41:39.760709 image/svg+xml Matplotlib v3.3.1, https://matplotlib.org/
2021-04-22T09:41:42.778239 image/svg+xml Matplotlib v3.3.1, https://matplotlib.org/
2021-04-22T09:41:46.119879 image/svg+xml Matplotlib v3.3.1, https://matplotlib.org/
In [26]:
data.lmPlot(dataType = 'matE', Erange = Erange, Etype = Etype, thres = 0.1, logFlag = False, selDims = {'Type':'L', 'it':1}, fillna = True, cmap = "vlag", pType = "phase", figsize = (15,8))   

# cmap = "cubehelix" )  # Quite good for showing resonance features.
Plotting data Xe_wf.orb19_4d_E1.0_15.0_286.0eV.inp.out, pType=phase, thres=0.1, with Seaborn
Plotting data XeF2_highRes_wf.orb21_A1G_E1.0_10.0_141.0eV.inp.out, pType=phase, thres=0.1, with Seaborn
Plotting data XeF2_highRes_wf.orb22_E1G_E1.0_10.0_141.0eV.inp.out, pType=phase, thres=0.1, with Seaborn
Plotting data XeF2_highRes_wf.orb24_E2G_E1.0_10.0_141.0eV.inp.out, pType=phase, thres=0.1, with Seaborn
2021-04-22T09:41:52.032513 image/svg+xml Matplotlib v3.3.1, https://matplotlib.org/
2021-04-22T09:41:55.040189 image/svg+xml Matplotlib v3.3.1, https://matplotlib.org/
2021-04-22T09:41:58.122883 image/svg+xml Matplotlib v3.3.1, https://matplotlib.org/
2021-04-22T09:42:01.431833 image/svg+xml Matplotlib v3.3.1, https://matplotlib.org/
In [ ]:
 

Versions

In [27]:
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter'])
Out[27]:
Thu Apr 22 09:42:03 2021 Eastern Daylight Time
OS Windows CPU(s) 32 Machine AMD64
Architecture 64bit RAM 63.9 GB Environment Jupyter
Python 3.7.3 (default, Apr 24 2019, 15:29:51) [MSC v.1915 64 bit (AMD64)]
epsproc 1.3.0-dev xarray 0.15.0 jupyter Version unknown
numpy 1.19.2 scipy 1.3.0 IPython 7.12.0
matplotlib 3.3.1 scooby 0.5.6
Intel(R) Math Kernel Library Version 2020.0.0 Product Build 20191125 for Intel(R) 64 architecture applications